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Abstract 

In  this  paper  we  give  stability  analysis  and  error  estimates  for  the  recently  introduced 
central  discontinuous  Galerkin  method  when  applied  to  linear  hyperbolic  equations.  A  com¬ 
parison  between  the  central  discontinuous  Galerkin  method  and  the  regular  discontinuous 
Galerkin  method  in  this  context  is  also  made.  Numerical  experiments  are  provided  to  vali¬ 
date  the  quantitative  conclusions  from  the  analysis. 
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1  Introduction 


In  this  paper  we  give  stability  analysis  and  provide  error  estimates  for  the  recently  introduced 
central  discontinuous  Galerkin  (DG)  method  [8]  for  linear  hyperbolic  equations,  and  compare 
the  central  DG  method  with  the  regular  DG  method  in  this  context. 

The  central  DG  method  is  designed  based  on  two  existing  techniques:  the  central  scheme 
framework  (e.g.  [9,  7])  and  the  DG  framework  (e.g.  [2,  4]).  It  uses  overlapping  cells  and  hence 
duplicative  information,  but  avoids  numerical  fluxes  (Riemann  solvers)  which  is  a  distinct 
advantage  of  central  schemes.  The  central  DG  scheme  also  avoids  the  excessive  numerical 
dissipation  for  small  time  steps,  common  to  some  of  the  earlier  central  schemes,  by  a  suitable 
choice  of  the  numerical  dissipation.  Being  a  variant  of  the  DG  method,  it  shares  many  of 
its  advantages,  such  as  compact  stencil,  easy  parallel  implementation,  etc.  The  central  DG 
method  performs  well  in  numerical  simulations  of  linear  and  nonlinear  scalar  and  systems  of 
conservation  laws  [8]. 

We  now  give  a  description  of  the  regular  DG  method  and  the  central  DG  method.  For 
simplicity,  we  consider  a  scalar  one  dimensional  conservation  law 

ut  +  f(u)x  =  0,  (x,  t)  G  [a,  b\  x  [0,  T]  (1.1) 

with  periodic  or  compactly  supported  boundary  conditions.  The  DG  and  central  DG  meth¬ 
ods  can  be  defined  for  nonlinear,  multi-dimensional  and  system  cases,  and  with  other  bound¬ 
ary  conditions  as  well,  see  [4,  8]. 

Let  {xj}  be  a  partition  of  [a,  b]  with  hj+ 1  =  xJ+\  —  Xj  and  h  =  max,,-  hJ+i.  The  mesh  is 
regular,  in  the  sense  that  max,,-  h-+i/  miry  h]+i  is  upper-bounded  by  a  fixed  constant  during 
mesh  refinements.  Denote  Xj+ 1  =  \{xj+ 1  +Xj),  Ij  =  (xj_  i ,  xJ+i ),  and  I-+ 1  =  (xj,Xj+ 1).  14 
is  the  set  of  piecewise  polynomials  of  degree  k  over  the  subintervals  {Ij}  with  no  continuity 
assumed  across  the  subinterval  boundaries.  Likewise,  ILy,  is  the  set  of  piecewise  polynomials 
of  degree  k  over  the  subintervals  {IJ+i}  with  no  continuity  assumed  across  the  subinterval 
boundaries. 
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The  regular  DG  method  is  defined  using  the  space  V/,  only.  The  semi-discrete  version  is 


as  follows.  Find  Uh(-,t)  G  14  such  that  for  any  (ph  G  14  and  for  all  j, 


dtuh(phdx  =  I  f{uh)dx(phdx-f[uh(x1,t),uh(x't1)t))(ph(x  i) 

J'  2 


j+r 


J+r 


+f  [Uh(x  ^t),uh(xj_i,t))  <ph(xj_i) 


(1.2) 


where  f(u~,u+)  is  a  monotone  numerical  flux,  namely  it  is  increasing  in  the  first  argument 
and  decreasing  in  the  second  argument,  or  symbolically  /('[',  J,);  it  is  consistent  with  the 
physical  flux  f(u,u )  =  f(u)\  and  it  is  at  least  Lipschitz  continuous  with  both  arguments. 
For  monotone  fluxes  suitable  for  the  DG  method,  see,  e.g.  [2],  For  systems  the  monotone 
numerical  flux  is  replaced  by  a  numerical  flux  obtained  through  an  exact  or  approximate 
Riemann  solver,  see,  e.g.  [10]. 

The  central  DG  method  is  defined  on  overlapping  cells  and  uses  both  spaces  14  and  114 • 
We  start  with  the  description  of  the  method  for  the  fully  discrete  version,  with  forward  Euler 
time  discretization  from  tn  to  tn+ i  =  tn  +  r  (here  the  time  step  r  could  change  with  n,  but 
for  simplicity  of  notations  we  will  use  a  constant  r).  The  scheme  is  defined  by  the  following 
procedure:  Find  uj]+1  G  14  and  v%+1  G  114,  such  that  for  any  iph  G  14  and  ijjh  G  114, 


uh+1[Phdx  =  0  I  v^iphdx  +  (1  -  9)  I  u^iphdx 


+T[  f  f(vh)dx<Phdx  -  f  (u£(xi+i))  <Ph{xj+i)  +  /  (<(4,-1))  <Ph(x^_  1)  )  (1-3) 


where  9  =  -1—  and  rmax  is  an  upper  bound  for  the  time  step  size  due  to  the  CFL  restriction, 
that  is,  rmax  —  ch  with  a  given  constant  CFL  number  c  dictated  by  stability.  Computa¬ 
tionally,  the  scheme  (1.3)-(1.4)  is  used  with  the  forward  Euler  replaced  by  a  Runge-Kutta 
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method  of  suitable  temporal  accuracy,  e.g.  the  SSP  Runge-Kutta  methods  in  [11,  5].  The 
time  step  r  can  be  chosen  arbitrarily  subject  to  the  stability  restriction  r  <  Tmax,  hence 
0  <  d  <  1.  If  0  =  1,  the  scheme  is  similar  in  spirit  to  the  original  central  scheme  [9]. 

We  now  take  the  limit  r  — >  0  to  obtain  the  semi-discrete  version  of  the  central  DG 
scheme:  Find  Uh(-,t )  G  14  and  Vh(-,t)  G  R4,  such  that  for  any  (fh  G  14  and  iph  G  114, 


dt uh  iphdx  = 


T, 


(vh  ~  uh)iphdx  +  /  f{yh)dxiphdx 


max  J  I, 
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~f  <Ph(xj+ 1)  +  /  (vh(Xj_  I,t))  <Ph(x+_  i)  (1.5) 


dtvh  i>hdx  = 


j+i 


tnax  J I  ,  i 
J+2 


(uh-vh)i>hdx+  f  (uh)dx,iphdx 


Ti+h 


-f  (uh(Xj+l ,t))  iph{xj+ 1)  +  /  («/.(%,«)) 


(1.6) 


Notice  that,  unlike  the  regular  DG  scheme  (1.2),  the  central  DG  scheme  (1.5)-(1.6)  does  not 
need  a  numerical  flux  to  define  the  interface  values  of  the  solution,  since  the  evaluation  of  the 
solution  at  the  interface  is  in  the  middle  of  the  staggered  mesh,  hence  in  the  continuous  region 
of  the  solution.  The  first  term  on  the  right  side  of  (1.5)  or  (1.6)  is  a  numerical  dissipation 
term.  This  will  become  clear  when  we  discuss  the  stability  of  the  scheme  in  Section  2.1.  In 
all  the  DG  schemes  (1.2),  (1.3)-(1.4)  and  (1.5)-(1.6),  the  initial  condition  is  taken  as  the  L 2 
projection  of  the  PDE  initial  condition  into  the  relevant  finite  element  space. 

The  organization  of  the  paper  is  as  follows.  In  Section  2,  we  first  analyze  the  L 2  stability 
and  give  an  a  priori  error  estimate  for  the  semi-discrete  central  DG  scheme  (1.5)-(1.6)  for 
the  linear  hyperbolic  equation,  using  similar  techniques  of  stability  and  error  analysis  for 
standard  DG  methods  [6,  3].  We  then  perform  a  Fourier  analysis  for  the  semi-discrete 
central  DG  scheme  (1.5)-(1.6)  for  the  linear  hyperbolic  equation  with  uniform  meshes  for 
piecewise  constant  and  linear  elements,  using  the  techniques  in  [13].  This  analysis  is  more 
explicit  and  allows  us  to  compare  the  errors  quantitatively  between  the  central  and  standard 
DG  schemes,  which  is  performed  in  Section  3.  In  Section  3  we  also  perform  numerical 
experiments  to  verify  such  quantitative  conclusions.  We  give  a  few  concluding  remarks  in 
Section  4. 
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2  Analysis  of  the  central  DG  scheme 


For  the  purpose  of  analysis  we  will  consider  the  linear  hyperbolic  equation,  namely  (1.1) 
with  f(u)  =  au  for  a  constant  a.  Our  analysis  can  be  easily  generalized  to  multi-dimensional 
linear  equations  and  hyperbolic  linear  systems. 

Without  loss  of  generality  we  consider  the  following  linear  hyperbolic  equation 


ut  +  ux  =  0,  (x,  t )  G  [a,  b]  x  [0,  T] 


(2.1) 


We  study  the  L2  stability  of  the  central  DG  scheme  (1.5)-(1.6)  for  the  equation  (2.1)  in 
Section  2.1,  and  compare  the  result  with  that  for  the  regular  DG  scheme  in  [6].  In  Section 
2.2  we  provide  an  L 2  a  priori  error  estimate  for  smooth  solutions,  and  compare  the  result 
with  that  for  the  regular  DG  scheme  in  [3].  In  Section  2.3  we  give  an  quantitative  error 
estimate  for  the  central  DG  scheme  for  polynomial  degree  up  to  1  using  Fourier  analysis, 
similar  to  the  technique  used  in  [12,  13]. 

2.1  L2  stability 


Theorem  2.1.  The  numerical  solution  Uh  and  Vh  of  the  central  DG  scheme  (1.5)-(1.6)  for 
the  equation  (2.1)  satisfies  the  following  L 2  stability  condition 


\jt  /  (K)2  +  (vh)2)dx  = 


Trr 


(Uh  -  Vhj2dx  <  0 


(2.2) 


Proof:  Taking  the  test  functions  (ph  =  Uh  and  'iph  =  Vh  in  (1.5)  and  (1.6)  respectively, 
summing  up  over  j,  observing  f{u)  =  u  and  the  periodic  (or  compactly  supported)  boundary 
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condition,  we  have 


1  d  rb 

2  dt  Ja 

1  fl 


1  max  Ja 


(K)2  +  (■ vh)2)dx 

(vhuh  -  ( uh )2  +  uhvh  -  ( vh )2)  dx  +  ^[J  vh  dxuhdx  +  J  uh  dxvhdx 


—vh  (xj+i,t)  uh  (xj+i,t)  +vh  (xj_i  ,tj  uh  (z+_i  ,tj 

~Uh  (Xj+i,t)  Vh  ( Xj+1,t )  +  Uh  ( Xj,t )  vh 

l  rh  rx3 

- /  (uh  ~  vhfdx  +  /  dx{uhVh)dx  + 


,  i 


dx(uhvh)dx 


J  Jxi-h 


-Vh  (xj+  i,tj  Uh  (xj+i  ,tj  +vh  (xj_i,tSj  uh 

-Uh  (xj+i,t)  Vh  (xj+1,t)  +  Uh  ( Xj,t )  vh  ( xj,t )] 

1 


(Uh  -  Vh)  dx  <  0 


Remark  2.1.  The  proof  of  Theorem  2.1  is  similar  to  the  proof  of  the  cell  entropy  inequality 
for  the  regular  DG  method  in  [6].  However,  we  cannot  prove  a  similar  L 2  stability  result 
for  the  central  DG  scheme  when  applied  to  the  nonlinear  scalar  conservation  law  (1.1),  even 
though  the  proof  of  Theorem  2.1  can  be  easily  generalized  to  multi-dimensional  central  DG 
schemes  for  linear  equations.  This  is  in  contrary  to  the  cell  entropy  inequality  for  regular 
DG  schemes,  which  holds  for  arbitrary  nonlinear  scalar  conservation  laws  [6]. 

Remark  2.2.  Theorem  2.1  indicates  that  the  energy  dissipation  term  is  —  j^(uh  —  Vh)2dx, 
that  is,  it  is  directly  related  to  the  difference  of  the  two  duplicative  representations  Uh  and 
Vh  of  the  solution  in  overlapping  cells.  In  contrast,  for  the  regular  DG  method,  the  energy 
dissipation  term  is  directly  related  to  the  jumps  of  the  numerical  solution  at  cell  interfaces. 

2.2  L2  a  priori  error  estimate 

In  this  subsection  we  use  the  standard  DG  techniques  [3]  to  obtain  an  a  priori  L 2  error 
estimate  for  the  central  DG  scheme. 
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Theorem  2.2.  The  numerical  solution  Uh  and  Vh  of  the  central  DG  scheme  (1.5)-(1.6)  for 
the  equation  (2.1)  with  a  smooth  initial  condition  u(-,0)  G  Hk+1  satisfies  the  following  L2 
error  estimate 

|| u  —  Uh\\2  +  || u  —  Vh\\2  <  Ch2k  (2.3) 

where  u  is  the  exact  solution  of  (2.1),  k  is  the  polynomial  degree  in  the  finite  element  spaces 
14  and  Wh,  and  the  constant  C  depends  on  the  (k  +  l)-th  order  Sobolev  norm  of  the  initial 
condition  ||m(-,0)||#m-i  as  well  as  on  the  final  time  t,  but  is  independent  of  the  mesh  size  h. 

Proof:  Let  us  first  introduce  the  standard  notation 


Bj(uh,  vh]  <ph,  iph)  =  /  dtuh  iphdx 


(vh  -  uh)(fhdx  -  /  vhdxtphdx 

+vh(xi+i,t)(ph(x~4_l)  -  vh{xj_i,t)(ph{x+_l)  (2.4) 

f  uhdx^hdx 


1 


d tvh  4 hdx 


3  +  i 


^max  J I  i 
J+rZ 


(uh  -  vh)^hdx 


j+i 


+uh{xj+1,t)^h(xj+l)  -  uh(xj.  I)vh(xj 


Clearly,  we  have: 


(2.5) 


for  all  j  and  all  (fh  G  14  and  iph  G  H4 •  It  is  also  clear  that  the  exact  solution  u  of  the  PDE 
(2.1)  satisfies 

Bj(u,u;tph,^h)  =  0  (2.6) 

for  all  j  and  all  <fh  £  14  and  4/i  £  H4-  Subtracting  (2.5)  from  (2.6),  we  obtain  the  error 
equation 

Bj(u-uh,u-vh;(ph,iph)  =  0  (2.7) 


for  all  j  and  all  (pt  G  14  and  iph  £  1T'4  - 

We  now  define  P  and  Q  as  the  standard  L 2  projection  into  14  arid  Wh  respectively.  That 
is,  for  each  j, 


(. Pw(x )  —  w(x))iph(x)dx  =  0  \/ifh  £  Pfc(/j) 


(2.8) 
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and 


(2.9) 


/  (Qw(x)  -  w(x))iph(x)dx  =  0  \/i>h  e  P k{Ij+±) 

Ji.M  i  2 

■?+!Z 

where  Pfc(/j)  and  Pfc(/J+i)  denote  the  spaces  of  polynomials  of  the  degree  up  to  k  in  the  cell 
Ij  and  the  cell  IJ+i  respectively.  Standard  approximation  theory  [1]  implies,  for  a  smooth 
function  w, 

\\(Pw(x)  —  w(x))\\  +  h1^2\\(Pw(x)  —  w(x))\\r ,  x  <  Chk+1  (2-10) 

■1  +  2 

and 

| \(Qw(x)  -  w(x))||  +  h1/2\\(Qw(x)  -  u'(x))||rJ.  <  Chk+1  (2.11) 

where  Tj  and  rj+i  denote  the  set  of  boundary  points  of  all  elements  Ij  and  Ij+1/2  respectively, 
and  the  positive  constant  4,  here  and  below,  solely  depending  on  w(x)  and  its  derivatives, 
is  independent  of  h. 

We  also  recall  that  [1],  for  any  Wh  G  14  or  Wh  G  114,  there  exists  a  positive  constant  C 
independent  of  wy,  and  h,  such  that 


\\dxwh\\  <  Ch  1\\wh\\]  ||wfe||r  <  Ch  1/2| 


wh\ 


(2.12) 


where  Y  =  Tj  or  I ]+i  ■ 
We  now  take: 


<Ph  =  Pu~  uh,  i\) h  =  Qu  -  vh 


(2.13) 


in  the  error  equation  (2.7),  and  denote 


ipe  =  Pu  —  u,  i\)e  =  Qu  —  u 


(2.14) 


to  obtain 


Bj((ph,  iph]  Wh,  'fph)  =  Bj(ipe,ipe]iph,iph). 


(2.15) 


For  the  left-hand  side  of  (2.15),  we  use  Theorem  2.1  to  conclude 


1  d 


^Bjiipk^hWh^h)  =  --j  /  ((<Ph)2  +  (iph)2)dx  + -  /  ((fh~^h)2dx.  (2.16) 

Z  CLT  J  n  ^max  J  a 


We  then  write  the  right-hand  side  of  (2.15)  as  a  sum  of  three  terms 


Bj{<pe,  ipe]  <ph,  iph)  =  B)  +  B2  +  B3 


3  '  3  1  3 


(2.17) 


where 


B)  =  —  f  (<fe  -  r)Vhdx  +  — 


(ip e  -  (fe)i>hdx  +  /  dt(peVhdx  +  /  dtipeiphdx 


1  max  J I .  i 


Bj=—  ipedx (phdx  -  /  Lpedxiphdx 


Bj  =  fpe(xj+i  ,t)(ph(x  i  ,t)-ipe(xj_i,t)(ph(x+  i ,  t)+^e(zj+i,  i)^/»(a;j+1,  *)-¥>e(zj,  t) 


and  we  will  estimate  each  term  separately. 
By  using  the  simple  inequality 


a/5  <  -(a2  +  (32), 


(2.18) 


the  L2  projection  property  (2.10)  for  (fe,  % ';e,  dttpe  and  d tipe,  and  the  fact  that  rmax  =  0(h), 


we  have: 


B 1  <  /  (yhfdx  +  /  (iph)2dx  +  C/r2fc+1 


(2.19) 


Likewise,  by  using  the  simple  inequality  (2.18),  the  L2  projection  property  (2.10)  for  tpe 
and  ipe,  and  the  first  inequality  in  (2.12)  for  and  iph,  we  have: 


B2  <  /  (<fh)2dx  +  /  (iphfdx  +  Ch 2 


(2.20) 


Finally,  by  using  the  simple  inequality  (2.18),  the  L 2  projection  property  (2.10)  for  <y3e 
and  ipe,  and  the  second  inequality  in  (2.12)  for  (ph  and  iph,  we  have: 


Bj  <  /  {(fh)2dx  +  /  (iph)2 dx  +  C7r2 


(2.21) 


Summing  up  (2.19),  (2.20)  and  (2.21)  over  j  and  combining  with  (2.16),  we  obtain  from 


(2.15) 


\jt  /  (M2  +  M2)dx  <  C  I  ((ph)2  +  (iph)2)dx  +  Ch2k 
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This,  together  with  the  approximation  result  (2.10),  implies  the  desired  error  estimate  (2.3). 


Remark  2.3.  The  error  estimate  of  Theorem  2.2  is  sub-optimal.  Numerically  we  observe 
the  optimal  (k  +  l)-th  order  accuracy,  see  [8]  and  also  the  numerical  results  in  next  section. 
In  contrast,  the  error  estimate  for  the  regular  DG  method  in  such  one  dimensional  and  also  in 
multi-dimensional  tensor  product  cases  is  optimal  [3].  The  major  technical  difficulty  leading 
to  this  loss  of  optimality  in  the  proof  is  the  B?  term.  For  the  regular  DG  method  this  term 
is  zero  for  both  the  regular  L 2  projection  and  for  a  special  projection  which  is  orthogonal  to 
polynomials  of  one  degree  lower  and  can  render  the  boundary  terms  in  B ?  also  to  vanish. 
However,  for  the  central  DG  method,  since  it  involves  i/)e  and  dx<fh  and  they  are  defined  by 
polynomials  in  different  cells,  it  is  impossible  to  make  the  Bj  term  vanish  no  matter  how  the 
projection  is  chosen,  although  a  special  projection  like  that  used  for  regular  DG  methods 
can  render  Bj  to  be  zero. 

Remark  2.4.  The  error  estimate  of  Theorem  2.2  can  be  easily  generalized  to  one  dimensional 
linear  hyperbolic  systems,  multi- dimensional  scalar  linear  hyperbolic  equations,  and  multi¬ 
dimensional  symmetric  linear  systems. 

2.3  A  quantitative  error  estimate  via  Fourier  analysis 

In  this  subsection  we  perform  a  Fourier  analysis  for  the  semi-discrete  central  DG  scheme 
(1.5)-(1.6)  for  the  linear  hyperbolic  equation  with  uniform  meshes  for  piecewise  constant 
and  linear  elements,  using  the  techniques  in  [13].  This  analysis  is  more  explicit  and  allows 
us  to  compare  the  errors  quantitatively  between  the  central  and  standard  DG  schemes. 

For  this  purpose  we  rewrite  the  scheme  (1.5)-(1.6)  for  the  linear  equation  (2.1)  as  a  finite 
difference  scheme  on  a  uniform  mesh.  Towards  this  goal  we  choose  the  degrees  of  freedom 
for  the  k-th  degree  polynomial  inside  the  cell  Ij  and  7J+i  respectively  as  the  point  values  of 
the  solution,  denoted  by 

u , ,  2 i-k  (t),  i  =  0 . k , 

J+2(fc+l)  V  ' 
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and 


v- ,  2i+i~k  (t),  i  =  0 . k , 

2(fc+l) 

at  the  k+1  equally  spaced  points 

/  2  i  —  k  \ 

d  +  2 (FTTyJ'1’  !  =  0""’A: 

and 

/.  2i  +  1  , 

d  +  ’  s  =  0-"’fc' 

The  schemes  written  in  terms  of  these  degrees  of  freedom  become  finite  difference  schemes 
on  a  globally  uniform  mesh  (with  a  mesh  size  h/{k  +  1)),  however  they  are  not  standard 
finite  difference  schemes  because  each  point  in  the  group  of  k+1  points  belonging  to  the  cell 
Ij  or  Ij+ 1  obeys  a  different  form  of  the  finite  difference  scheme. 

To  be  more  specific,  we  first  concentrate  on  the  simplest  piecewise  constant  k  =  0  case. 
For  this  case,  we  choose  the  degrees  of  freedom  as  the  point  values  at  the  N  uniformly  spaced 
points 

Uj(t),  j  =  0,  ...,7V-  1. 
or 

vj+ §(*),  j  =  0,  ...,7V-  1. 

The  solution  inside  the  cell  Ij  or  IJ+ 1  is  then  represented  by 

uh(x,t)  =  Uj(t)<p°h(x) 

or 

Vh(x,t)  =  vj+i(t)ip%(x) 

where  +Gh[x)  is  the  constant  function  which  equals  1  inside  Ij,  and  similarly  x )  is  the 
constant  function  which  equals  1  inside  7?+i .  With  this  representation,  taking  the  test 
functions  ifh  as  ip°h,  and  iph  as  -0°,  respectively,  we  obtain  easily  the  finite  difference  schemes 
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corresponding  to  the  central  DG  scheme: 


U3  = 


^LUj  +  (db  + 1) v ^  +  (zb  - 1) 


lU  ~  ■dil'^  +  (s~  +  s)“j  +  (s~_s)"J+1  (2-22) 


for  j  =  0, N  —  1.  Here  v!  and  v'  denote  the  time  derivatives  of  u  and  v.  The  scheme  can 
be  rewritten  into  a  more  compact  form 


Uj  )+c  Uj+1 


(2.23) 


0  -^  +  4 

^  1  max  u> 

0  0 


i  i  _  i 

Tmax  2Tmax  h 

_  +  I _ i_ 

r  h  ^ 


0  0 
0 

ix  h 


(2.24) 


We  now  perform  the  following  standard  Fourier  analysis  for  the  finite  difference  scheme 
(2.23)-(2.24).  This  analysis  depends  heavily  on  the  assumption  of  uniform  mesh  sizes  and 
periodic  boundary  conditions.  We  make  an  ansatz  of  the  form 


Uj(t) 

V3+h(t) 


Um(t ) 

Vm(t) 


(2.25) 


and  substitute  this  into  the  scheme  (2.23)-(2.24)  to  hnd  the  evolution  equation  for  the  coef¬ 


ficient  vector  as 


Mm(t)  \  Q,  /  x  f  um(t ) 


(2.26) 


where  the  amplihcation  matrix  G(m,  h )  is  given  by 


G{m,h)  =  Ae~imh+  B  +  C  e* 


(2.27) 


with  the  matrices  A,  B  and  C  defined  by  (2.24).  The  two  eigenvalues  of  the  amplihcation 
matrix  G(m,  h )  are 


X  ■  mh 

- ae  2  ; 


1  i  m/t 

- V  ae  2 


(2.28) 
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where  a  =  ^ - b  j-)e  imh  +  75—^ - f-  The  general  solution  of  the  ODE  (2.26)  is  given  by 

v  -^Tmax  7  ^Tmax  v  7 

{tm)  =  a^+a^  <2'29) 

where  the  eigenvalues  Ai  and  A2  are  given  by  (2.28),  and  Vj  and  V2  are  the  corresponding 
eigenvectors  given  by 


Vj  = 


V2  = 


with  £  =  mh.  For  accuracy  we  look  at  the  low  modes,  in  particular  at  m 
given  initial  condition 

«;(0)  =  e“',  vj+l(0)=ex^ 


(2.30) 
1.  To  fit  the 


(2.31) 


whose  imaginary  part  is  the  initial  condition 


Wj(0)  =  sin(xi),  vj+i  (0)  =  sin(xj+i), 

we  require,  at  t  —  0, 


hence  we  obtain  the  coefficients  a\  and  a2  in  (2.29)  as 

„  i  k 

a\  —  0,  a2  =  e  2. 


(2.32) 


(2.33) 


We  remark  that  the  usual  way  of  taking  initial  conditions  in  a  finite  element  method  is  via 
an  L 2  projection,  not  by  a  point  value  collocation  (2.32),  however  we  have  verified  that  this 
does  not  affect  the  final  results  in  the  analysis  in  this  paper.  We  thus  have  the  explicit 
solutions  of  the  scheme  (2.23)-(2.24)  with  the  initial  condition  (2.31),  for  example 

Uj(t)  =  a2eiXj+X2t-^  (2.34) 

with  the  eigenvalue  A2  given  by  (2.28)  with  m  —  1  and  the  coefficient  a2  given  by  (2.33).  By 
a  simple  Taylor  expansion,  we  obtain  the  imaginary  part  of  Uj(t )  to  be 

Im{uj(t)}  =  sin(xj  —  t)  —  f  - sin (xj  —  t)  \  h  +  0(h2) 


13 


(2.35) 


where  c  = 


=  0(1)  is  the  maximum  CFL  number.  This  is  clearly  consistent  with  the 


exact  solution  to  first  order  accuracy.  We  can  similarly  check  the  first  order  accuracy  of 
Im{vj+i(t)}. 

We  now  repeat  this  analysis  for  the  piecewise  linear  k  —  1  case.  The  solution  inside  the 


cell  Ij  or  I-+ 1  is  then  represented  by 


uh{x,t)  =  +  uj+i{t)(pl{£) 


or 


Vh{x,t)  =  Vj+1  +  vj+*  Vh(£) 

where  <?£(£)  =  -f  +  ip2h(£)  =  £  +  ip\{£)  =  -f  +  §  and  ipl(£)  =  £  -  with  £  = 

2{x^x^  ■  With  this  representation,  taking  the  test  functions  iph  and  ip h  as  ip\,  (f>\  and  ip\,  ip \ 
respectively,  and  inverting  the  small  4x4  mass  matrix  by  hand,  we  obtain  easily  the  finite 
difference  scheme  corresponding  to  the  central  DG  scheme  as 


/ 


V1 

J  4 

u 


\ 


i+i 

u 

\  v'i+l  J 


=  A 


(  ui-\  \ 

Ua_  3 

J  4 

V~_3 
J  4 

V  vi-\  ) 


(  uJ-\  ^ 


B 


u 


4+4 

'•i+i 

V  / 


+  o 


\ 


(  mj+! 
mj+! 
Vt+I 


V  V4+f  / 


(2.36) 


with 


A  = 


5  = 


C  = 


(  0  0 
0  0 
0  0 
V  o  o 


1  +  A 

16Tmax  4ft 


16th 


4ft 


13 


16r, 


- L  \ 

4ft  \ 

_ +  A 

16Tmax  4  h 

0 


3“ 


-1 

Tmax 

0 

1  +  A 

16rmax  4  ft 


-1 


16th 


0 

0 


16^a 

16Tma 


13 


0 

-1 

Tmax 


4ft 


16rmax  4ft 
3  I  _5_ 
16rmax  ‘  4ft 


0 

0 


__  _5_ 
4  ft 

+  A 

~  4ft 


-1 


16r, 


_jna 


16th 


16Tmax 

13 

16Tmax 

-1 


+ 


_5_ 

4  ft 

_j_ 

4ft 


Tmax 

0 


0  0  \ 
0  0 
+  rh  o  o 
4  o  o/ 


-i 


16Tmax 
1  _ 

16Tmax 

0 

-l 


+  ik  \ 


_5_ 

4ft 


(2.37) 
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We  make  an  ansatz  of  the  form 


( Vl(*)  \ 

(  V-i(^)  \ 

VI  (*) 

V|W 

VI  (*) 

V|W 

V|W  / 

V 

(2.38) 


and  substitute  this  into  the  scheme  (2.36)-(2.37)  to  find  the  evolution  equation  for  the  coef¬ 
ficient  vector  as 


/  “L-iW  \ 


u'  i(f) 
m,+|v  > 


\ 


( 


=  G(k,  h ) 


u. 


^,|W 

V'm, 

where  the  amplification  matrix  G(m,  h )  is  given  by 


/ 


,-iW  t 
,++) 
VjW 

V  v§W  ) 


(2.39) 


G(m,/i)  =  Ae~imh+  B  +  Ceimh. 


(2.40) 


with  the  matrices  A,  B  and  C  defined  by  (2.37).  The  eigenvalues  of  the  amplification  matrix 
G(m,  h )  are 

Ai  =  — — (— 1  +  ~^-e~imh \J —a\  —  a2) 

Vax  * 

A2  =  - (— 1 - —e  imh\/ —a.\  —  a2) 

7~max  * 

A3  =  — (-l  +  yWV+^) 

^max  * 

A4  =  — — (  — 1  —  ye'^vW) 

^max  O 

where 

cq  =  (1  —  4c  —  52c2)eimh  —  2(11  —  52c2)e2imh  +  (1  +  4c  —  52c2)e3imh, 
a2  =  eimh(l  +  eimh  +  10c(l  -  eimh ))  [-3  -  12c  +  4c2  (2.41) 

+2(21  -  4 c2)eimh  +  4(— 3  +  12c  +  4c2)e2*m/l]1/2 


with  c  =  being  the  maximum  CFL  number. 
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The  general  solution  of  the  ODE  (2.39)  is  given  by 
/  um  -i(t)  \ 

=  a i  eXlt  Ei  +  a2  eA2‘  E2  +  a3  eAst  E3  +  a4  eA4'  E4,  (2.42) 

V  V§(*)  / 

where  the  eigenvalues  A4  ,  A2,  A3  and  A4  are  given  by  (2.41),  and  V\  ,  E2,  E3  and  E4  are  the 


corresponding  eigenvectors  given  by 


Ei  = 


V,  = 


E3  = 


E4  = 


/ 


s/2e 


V~Ql^a2(«7+2Q2) 


05(06— 02) 


\ 


y/-2(a1+a2)(l+eimh-10c(-l+eimh)) 

QlQ — 0=2 
OL4.CX.2  ~I~Q=3 
£*5(0=6—  CK2) 

1 

^/2e~lmhy/— 01—02(07+202) 


) 


_ 015(06— 02) 

y/-2(ai+a2)(l+eimh-10c(-l+eimh)) 

OQ—Q2 

0402+0.3 

05(06—02) 

1 


+2e 


^\/— 01+02+07+202) 


05(06+02) 


^2(-ai+a2)(l+eirnh-Wc(-l+eimh)) 

O6+O2 
—  Q4Q2+P3 
05(06+02) 

1 

_ \f2e~lvllh  \J — ai  +02  (—07+202) 

05(06+02) 


\ 


yj  2(—ai+O2)(l+eirnh  —  10c(—l+eirnh)) 

06+02 
— 0402+03 

05(06+02) 

1 


(2.43) 


where  c  =  — ^  is  still  the  maximum  CFL  number,  <+1  and  a2  are  given  by  (2.41),  and  the 


remaining  a's  are  given  by 


a3  =  -(7  +  104c  +  332c2  -80c3)eimh  +  (77  +  896c  +  532c2  -  240c3)e2im/l 
+(79  -  856c  -  68c2  +  240c3)e3im/l  -  (5  -  64c  +  132c2  +  80 c3)e4im/\ 
a4  =  -1  -20c+  (-3  +  20c)eim/\  (2.44) 

a5  =  -1  -4c+  (13  +  4c)eim/l, 

a6  =  (1  +  10  c)2eimh  +  (2  -  200c2)e2im/l  +  (1  -  10c)2e3im\ 
a7  =  (1  +  10 c)eimh  -  20 c2e2imh  -  (1  -  10 c)2e3imh. 
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We  again  look  at  the  low  modes  to  determine  accuracy.  In  particular  we  look  at  m  —  1. 


To  fit  the  given  initial  condition 


IX.  1 

Uj_i( 0)  =  e 


u 


3+ 


i(0)  = 


IX  .  .  1 


”j+j 


i(0)  = 


IX  .  ,  1 


D'+ 


3(0)  = 


3 


(2.45) 


whose  imaginary  part  is  our  initial  condition  for  (2.1),  we  require,  at  t  —  0, 


(  \ 

e*  <r 

V  i(°) 

i  - 

e  4 

v§(°)  ) 

•  3  h 

\  e  4  / 

This  gives  us  the  coefficients  a4,  a2 ,  03  and  0.4  in  the  solution  (2.42).  We  thus  have  the 
explicit  solution  of  the  scheme  (2.36)-(2.37)  with  the  initial  condition  (2.45),  for  example 

u,  1  =  aie^+Ai*-^I/i  +  a2eixi+X2t~i!tV2  +  a3eix*+x  3t~^V3  +  a4eiX3+Mt~^  V4  (2.46) 

J  4 


with  the  eigenvalues  Ai,  A2,  A3,  A4  given  by  (2.41)  and  the  eigenvectors  V\ ,  V2,  V3,  V4  given  by 
(2.43)  with  m  =  1,  and  the  coefficients  a4,  a2,  a3  and  a4  obtained  above  by  fitting  the  initial 
condition.  Through  a  simple  Taylor  expansion,  we  obtain  the  imaginary  part  of  u,  i(t)  to 

J  4  v 

be 

sin(  x  •  1  —  t ) 

Im{uj_i(t)}  =  sin(^_i  -  t) - g - h2  +  0(h3)  (2.47) 

for  a  fixed  choice  of  rmax  =  OAh.  Results  for  other  choices  of  rmax  also  indicate  the  same 
second  order  accuracy.  The  results  for  Im{Uj+i(t)},  Im{Vj+i(t)}  and  Im{v^+z{t)}  are 

similar. 

In  principle  this  analysis  can  be  performed  for  higher  order  polynomials  in  the  central 
DG  scheme,  however  the  algebra  becomes  prohibitively  complicated. 


3  A  comparison  between  the  central  DG  and  standard 
DG  methods 


In  [13],  results  to  similar  to  those  in  (2.35)  and  (2.47)  were  obtained  for  the  regular  DG 
scheme  (1.2)  applied  to  the  linear  equation  (2.1)  with  an  upwind  numerical  flux.  For  the 
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piecewise  constant  k  —  0  case,  the  result  for  the  regular  DG  scheme  is 

Im{uj(t)}  =  sin(xj  —  t)  —  sin(xj  —  t)^j  h  +  0(h2),  (3.1) 

and  for  the  piecewise  linear  k  =  1  case,  it  is 

sin(  x  •  i  —  t ) 

Im{uj_i(t)}  =  sin(^_i  -  t) - ^ - h 2  +  0(h3)  (3.2) 

We  are  now  in  a  position  to  make  a  quantitative  comparison  between  the  regular  and 
central  DG  schemes.  For  the  piecewise  constant  k  =  0  case  we  have  the  following  conclusions: 

1.  The  semi-discrete  versions  of  the  regular  and  central  DG  schemes  are  both  stable. 
When  discretized  with  the  first  order  forward  Euler  method,  the  CFL  numbers  for  the 
DG  method  and  for  the  central  DG  method  are  1.0  and  0.5,  respectively  (this  can  be 
veriEed  by  an  easy  von  Neumann  analysis).  When  discretized  with  the  second  order 
nonlinearly  stable  Runge-Kutta  method  [11],  the  CFL  numbers  for  the  DG  method 
and  for  the  central  DG  method  are  1.0  and  0.87,  respectively.  Thus  the  central  DG 
method  has  a  smaller  CFL  number. 

2.  They  are  both  first  order  accurate.  By  a  comparison  of  (2.35)  and  (3.1),  the  leading 

errors  for  the  central  and  regular  DG  methods  for  the  first  mode  (i.e.  for  the  sin(a;) 
initial  condition)  have  a  ratio  l/(4c).  That  is,  the  central  DG  method  has  a  smaller 
error  than  the  standard  DG  method  on  the  same  mesh  when  c  =  >  4. 

We  now  compute  the  DG  and  central  DG  solutions  to  (2.1)  with  u(a:,0)  =  sin(a;)  as  the 
initial  condition  and  with  periodic  boundary  conditions,  up  to  t  =  25  (about  four  periods 
later  in  time),  to  verify  the  quantitative  comparison  above.  In  our  computation  we  take 
rmax  =  0.8 h,  namely  c  =  0.8.  We  take  a  small  time  step  r  =  O.Olh  to  reduce  the  effect  from 
the  time  discretization.  In  order  to  be  consistent  with  the  error  analysis  above,  the  errors 
are  computed  for  Uh  at  the  points  Xj.  The  L 2  and  L°°  errors  and  order  of  accuracy  of  the 
central  DG  and  regular  DG  methods  are  listed  in  Tables  3.1  and  3.2  respectively.  We  also 
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Table  3.1:  L2  and  L°°  errors  and  orders  of  accuracy  of  the  first  order  central  DG  method. 


numerical  solutions 

predicted  by  analysis 

h 

order 

L°°  error 

order 

order 

L°°  error 

order 

2tt/80 

1.88E-01 

— 

2.65E-01 

— 

2.17E-01 

— 

3.07E-01 

— 

2tt/160 

1.01E-01 

0.90 

1.43E-01 

0.90 

1.09E-01 

1.00 

1.53E-01 

1.00 

2tt/320 

5.23E-02 

0.95 

7.40E-02 

0.95 

5.42E-02 

1.00 

7.67E-02 

1.00 

2vr/640 

2.67E-02 

0.97 

3.77E-02 

0.97 

2.71E-02 

1.00 

3.84E-02 

1.00 

2tt/1280 

1.35E-02 

0.99 

1.90E-02 

0.99 

1.36E-02 

1.00 

1.92E-02 

1.00 

Table  3.2:  L2  and  L°° 


errors  and  orders  of  accuracy  of  the  first  order  regular  DG  method. 


numerical  solution 

predicted  by  analysis 

h 

L1  2  error 

order 

L°°  error 

order 

L 2  error 

order 

L°°  error 

order 

2vr/80 

4.42E-01 

— 

6.25E-01 

— 

6.94E-01 

— 

9.81E-01 

— 

2tt/160 

2.74E-01 

0.69 

3.88E-01 

0.69 

3.47E-01 

1.00 

4.91E-01 

1.00 

2vr/320 

1.54E-01 

0.83 

2.18E-01 

0.83 

1.73E-01 

1.00 

2.45E-01 

1.00 

2vr/640 

8.17E-02 

0.91 

1.16E-01 

0.91 

8.68E-02 

1.00 

1.23E-01 

1.00 

2vr/1280 

4.21E-02 

0.96 

5.95E-02 

0.96 

4.34E-02 

1.00 

6.14E-02 

1.00 

list  the  predicted  errors  by  the  analysis,  namely  the  leading  terms  in  the  Taylor  expansions 
in  (2.35)  and  (3.1)  in  these  tables.  We  can  see  that  the  predicted  errors  and  the  actual  errors 
are  very  close,  validating  our  quantitative  analysis  in  Section  2.3. 

Likewise,  for  the  piecewise  linear  case  k  —  1,  we  have  the  following  conclusions: 

1.  The  semi-discrete  versions  of  the  regular  and  central  DG  schemes  are  both  stable. 
When  discretized  with  the  second  order  nonlinearly  stable  Runge-Kutta  method  [11], 
the  CFL  numbers  for  the  DG  method  and  for  the  central  DG  method  are  0.33  and 
0.45,  respectively  (this  can  be  verified  by  an  easy  von  Neumann  analysis).  Thus  the 
central  DG  method  has  a  larger  CFL  number. 

2.  They  are  both  second  order  accurate.  The  central  DG  method  has  different  leading 
errors  bh2  sin(xJ-_i  —  t)  for  different  c  =  with  the  sin(x)  initial  condition.  We  list 
the  corresponding  b  with  different  values  of  c  in  Table  3.3. 
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Table  3.3:  The  leading  error  term  bh2  sin(x-_i  —  t)  for  the  central  DG  method  for  different 

J  4 

p  _  Tmax 

(  ~  h 


c 

0.2 

0.3 

0.4 

0.5 

0.6 

0.8 

~7T 

l 

l 

l 

l 

l 

2 

u 

_ 

40 

30 

24 

_ 2 

_ 

Table  3.4:  L2  and  L°° 


errors  and  orders  of  accuracy  of  the  second  order  central  DG  method. 


numerical  solution 

predicted  by  analysis 

h 

L2  error 

order 

L°°  error 

order 

L 2  error 

order 

L°°  error 

order 

2tc/20 

1.20E-02 

— 

1.37E-02 

— 

4.10E-03 

— 

4.11E-03 

— 

2vr/40 

1.53E-03 

2.98 

1.93E-03 

2.82 

4.11E-04 

2.00 

4.11E-04 

2.00 

2vr/80 

1.91E-04 

3.00 

2.94E-04 

2.72 

1.03E-04 

2.00 

1.03E-04 

2.00 

2vr/160 

2.57E-05 

2.89 

4.96E-05 

2.57 

2.57E-05 

2.00 

2.57E-05 

2.00 

2vr/320 

6.43E-06 

2.00 

9.42E-06 

2.40 

6.43E-06 

2.00 

6.43E-05 

2.00 

By  a  comparison  of  Table  3.3  and  (3.2),  the  leading  errors  for  the  central  and  regular 
DG  methods  for  the  first  mode  (i.e.  for  the  sin(x)  initial  condition)  are  equal  when 
c  =  Ists^  =  0.5.  The  central  DG  method  has  a  smaller  error  than  the  standard  DG 

h 

method  on  the  same  mesh  when  c  <  0.5. 

We  now  compute  the  DG  and  central  DG  solutions  to  (2.1)  with  a  u(x,  0)  =  sin(x)  initial 
condition  and  periodic  boundary  conditions,  up  to  t  =  25  (about  four  periods  later  in  time), 
to  verify  the  quantitative  comparison  above.  In  our  computation  we  take  rmax  =  0.2 h,  namely 
c  =  0.2.  We  take  a  small  time  step  r  =  O.Ol/i  to  reduce  the  effect  from  the  time  discretization. 
In  order  to  be  consistent  with  the  error  analysis  above,  the  errors  are  computed  for  Uh  at  the 
points  Xj_  i  and  Xj+ 1.  The  L 2  and  L°°  errors  and  order  of  accuracy  of  the  central  DG  and 
regular  DG  methods  are  listed  in  Tables  3.4  and  3.5  respectively.  We  also  list  the  predicted 
errors  by  the  analysis,  namely  the  leading  terms  in  the  Taylor  expansions  in  Table  3.3  and 
(3.2)  in  these  tables.  We  can  see  again  that  the  predicted  errors  and  the  actual  errors  are 
very  close,  validating  our  quantitative  analysis  in  Section  2.3. 
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Table  3.5:  L2  and  L°°  errors  and  orders  of  accuracy  of  the  second  order  regular  DG  method. 


numerical  solution 

predicted  by  analysis 

h 

L2  error 

order 

L°°  error 

order 

L2  error 

order 

L°°  error 

order 

2tt/20 

1.06E-02 

— 

1.46E-02 

— 

4.10E-03 

— 

4.11E-03 

— 

2vr/40 

1.34E-03 

2.98 

2.36E-03 

2.63 

1.03E-03 

2.00 

1.03E-03 

2.00 

2vr/80 

2.57E-04 

2.39 

4.24E-04 

2.47 

2.57E-04 

2.00 

2.57E-04 

2.00 

2vr/160 

6.42E-05 

2.00 

8.53E-05 

2.32 

6.43E-05 

2.00 

6.43E-05 

2.00 

2vr/320 

1.61E-05 

2.00 

1.82E-05 

2.19 

1.61E-05 

2.00 

1.61E-05 

2.00 

4  Concluding  remarks 

We  have  performed  an  L 2  stability  analysis  and  an  a  priori  error  estimate  for  the  recently 
introduced  central  discontinuous  Galerkin  method  when  applied  to  linear  hyperbolic  equa¬ 
tions.  We  have  also  performed  a  Fourier  type  error  analysis  which  is  more  quantitative  and 
allows  us  to  make  a  comparison  with  the  regular  discontinuous  Galerkin  method.  It  is  verified 
that,  even  though  the  central  discontinuous  Galerkin  method  uses  duplicative  representa¬ 
tion  of  the  solution,  hence  involves  twice  the  computational  cost  and  storage  requirement 
than  the  regular  discontinuous  Galerkin  method,  it  is  more  accurate  for  certain  choices  of  a 
dissipation  parameter  for  the  same  mesh.  The  stability  analysis  and  error  estimates  do  not 
seem  to  be  easily  generalizable  to  nonlinear  hyperbolic  equations.  Further  analysis  in  this 
direction  is  needed.  A  comprehensive  comparison  of  the  numerical  performance  of  the  central 
discontinuous  Galerkin  method  and  the  regular  discontinuous  Galerkin  method  for  nonlinear 
multi-dimensional  systems  of  hyperbolic  conservation  laws  would  also  be  very  useful. 
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